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Chapter One 
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^ Wave propagation in one-dimension: 

> l\/lethods and applications to complex and 

;z; fractal structures 



O 



E. Akkermans*^^^ G. V. Dunne*^^^ and E. Levy^^) 

(^) Department of Physics, Technion, Israel Institute of Technology, 32000, Haifa, Israel 
(^) Department of Physics, University of Connecticut, Storrs CT 06269, USA 

O This chapter is a pedagogical review of methods and results for studying wave propagation 

C/^ in one-dimensional complex structures. We describe and compare the tight-binding, scat- 

O tering matrix, transfer matrix and Riccati formalisms. We present examples for transport 

(/^ through finite-sized layered dielectric systems with periodic, quasi-periodic, fractal, disor- 

dered, and random structure, illustrating how can spatial structure affect the spectrum of 

modes as well as the local mode intensity. 

<N 
> 

Qn 1.1. Introduction 

O 

Both classical and quantum waves are highly sensitive probes of the details of the 
physical and geometric structure of the medium in which they propagate. Peri- 

^ odic structures produce spectra with transmission bands, leading to the remark- 

able physics of photonic band gaps.^ Other kinds of structure, quasi-periodic,^'^ 
• • fractal, disordered or random,^"^ lead to other distinct spectral features, some of 

. ^ which suggest certain advantages for physical and engineering applications. Im- 

portant questions include the frequency distribution of modes, as well as their 
^ spatial localization. Recent experimental advances suggest that various tj^pes of 

complex media can now be fabricated with high precision. This review is designed 
to summarize some basic theoretical tools for designing layered structures in order 
to produce desired properties of the mode spectra and the local mode densities. 

We describe the density of states, the counting function (the integrated density 
of states), the transmission probability, and the local field amplitude and intensity, 
each of which has a direct physical meaning. We compare and contrast different 
computational methods: the tight-binding, scattering matrix, transfer matrix, and 
Riccati methods, and we present a variety of results to illustrate the effect of the 
layered structure type on the aforementioned physical quantities. 
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1.2. Wave equations 

The general aspects of coherent propagation are common to a wide variety of 
waves which propagate in scattering media. This notwithstanding, each type of 
wave exhibits its own characteristic behavior. In this section, we present several 
examples of wave equations, and we study two important classes in greater de- 
tail : the Helmholtz equation, which describes scalar wave propagation suited for 
electromagnetic TE or TM modes propagating in a dielectric, and the Schrodinger 
equation associated with a non interacting electron gas (weakly disordered metals 
or semiconductors). 

1.2.1. Helmholtz equation 

The case of electromagnetic waves is special, for several reasons. It is probably one 
of the earliest examples where changes in wave phase coherence due to passage 
through a random medium was examined. In the beginning of the twentieth cen- 
tury, very precise studies were carried out on electromagnetic wave propagation 
through diffusive media, specifically the atmosphere. From a conceptual view- 
point, this problem stimulated the community working in the theory of probabil- 
ity, who regarded it as a new field for the application of methods developed for 
the study of Brownian motion.^ For the atmosphere, the description in terms of a 
static disordered medium is not appropriate. For many other cases, however, the 
description in terms of static disorder described by a time-independent potential 
works well, and it is this case that we consider in the present section. 

For a TE mode propagating in the x direction along a medium with spatially 
varying dielectric function Maxwell's equations for the electromagnetic field 
amplitudes, Ez = ^p{x)e~^^^ and Hy = become two coupled equations:^ 

where /cq = noj/ c, namely for an average dielectric of refractive index n = V^/^o 
and average dielectric constant e (See for instance,^ Chap.2). This reduces to a 
scalar Helmholtz equation for 

-jS-^o^^^M^)=koH^) (1-2) 



where v{x) = — 1 = n^{x) — 1, with refractive index n{x). 



1 .2.2. Schrodinger equation 

The Schrodinger wave equation for a particle of mass mm a. one dimensional po- 
tential V{x) also has a Helmholtz-like form: 

This can be formally written a.s H^p = k^\p. 
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1.3. Tight binding formalism | 3 

The two examples we have discussed are not the only ones to exhibit effects 
related to coherent propagation in complex media. In fact, these effects are com- 
mon to all wave phenomena (quantum, optical, hydrodynamic, etc), independent 
of dispersion relation and space dimension, provided that there are no nonlinear 
effects. Indeed, these may hide complexity effects. Moreover, nonlinear equations 
often have special solutions (solitons, vortices,...) whose stability is ensured by a 
topological constraint which is very difficult to destabilize by means of a disor- 
dered potential. The competitive role of disorder and nonlinearity is important 
but is still relatively poorly understood.^^ 



1.3. Tight binding formalism 

For certain applications, a discrete version of the Helmholtz and Schrodinger equa- 
tions may be more appropriate or more tractable numerically. The tight-binding 
version of the Schroding er eq uation is widely used, for example in condensed mat- 
ter physics. ^'^^ We write \l3\ in the general form: 

Y^-vix)^p^-k^^p (1.4) 
Discretize space into a succession of identical layers of thickness A: 

^ [%+i{k) - 2xpn{k) + %-i{k)] + [k^ - vn] M^) = (1-5) 
In the limit A ^ 0, the previous equation can be rewritten as: 

^'P„+i{k) + ^%-i{k) - (l-k^A^) e^'^"Mk) (1.6) 

Defining the dimensionless quantities: e = 1 — fc^A^, (pn{e) = e^^^'^^'^xpn and 
= e^'^^^/^^'^^^+^^+i^ the wave equation takes the canonical discrete form: 

(^) + (pn-1 {e) =2e (pn (1.7) 

often referred to as non-diagonal tight-binding description. The equivalent diagonal 
tight-binding form is: 

(Pn+l (e) + (pn-1 {e) +Vn(pn=2e(pn (1.8) 

These two equivalent forms define three-term recursion relations,^^ whose solu- 
tions are families of polynomials. They show up in various problems such as the 
master equation, random walks, and laplacian problems. An important issue is 
to determine the corresponding family of properly normalized orthogonal poly- 
nomials and their spectral distribution. For example, the solutions for a periodic 
potential are in terms of Chebychev polynomials, as discussed below in Section 
[TO 

For a system of length (number of sites) N, the tight-binding equation can be 
written in a matrix form Hjs/O = e O. The system has thus N eigenenergies, we 
shall denote by ei, . . . , e^. The counting function, N'{e), or integrated density of states, 
is defined as the fraction of eigenenergies that are smaller than a given energy e: 

1 N 

^f{e)^^Yl^{e-ei) (1.9) 

1=1 
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where 6{x) is the Heaviside step function, equal to 1 (resp. 0) for x > (resp. 
X < 0). For large enough N, the counting function is independent of the choice 
of boundary conditions. The counting function is usually a well defined and con- 
tinuous function of energy Its derivative, p{e) = dJV/de, if it exists, is called the 
density of states (or density of modes). It is a sum of Dirac delta functions, 

pW = ]^E^(^-^0- (1-10) 

As we shall see, for fractal systems, the counting measure dM{e) exhibits a singular 
behavior related to divergences of the density of states. 

Remark : There is a connection between the counting function and the zeroes of 
the wave function in one-dimensional systems. Consider for simplicity the form 
( |1.8[ ) of the tight-binding wave equation for N sites. Assuming zero boundary 
conditions (Dirichlet), (p^ = (pN+i = 0/ and (pi = 1 as a normalization, we have, 
(p2 = e — vi, (p3 = {e — vi){e — V2) — 1, and so on, so that, 

N 

(pN+i=Yl('-ei). (1.11) 

i=l 

namely the zeroes of the wave function are related to the location of eigen-energies. 



The tight binding description allows to make interesting connections to other 
descriptions such as the transfer matrix approach (see sections [L5] and 1.8.3[ ). A 



simple illustration for the connection is the simple case of the laplacian evolution 
in the absence of scattering potential. It corresponds to ( |1.4| for v{x) = 0. For 
the tight binding problem it corresponds to taking tn,m = 1 for all sites, namely to 
solve the three term recursion: 

(Pn+i{e) =2e(Pn-(pn-i{e). (1.12) 

The solution of this specific recursion is well known and given by the Cheby- 
shev polynomials of the first kind. For —l<e<l, the solutions are (pn{e) = 
(pn{cosO) = Tn{cosO) = COS n^, where e = cos^, and G [—re, re] is the so called 
Bloch angle. Their spectral distribution is characterized by the density of states: 

which corresponds to a well known results both in tight-binding and for the in- 
variant spectral measure of Chebyshev polynomials.^^ 



1 .3.1 . Tight'binding approach to the Helmholtz equation 

The tight binding version of the scalar wave (Helmholtz) equation is obtained 
along the same lines. To derive it, consider instead of the continuous variation 
of the refractive index n(x), a succession of N dielectric layers extending from 
< X < xi, xi < X < X2/ • • • / ^N-i < ^ < ^N- For each layer of thickness there 
is a constant refractive index n/. We can then rewrite \i.6} as: 

^i+i = K^i - Mi-i (1-14) 
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1.4. Scattering matrix formalism 5 



where 



fii sin j6v 1 1 



(1.15) 



and j6f = koUili. 

1.4. Scattering matrix formalism 

Propagation through a one-dimensional structure can be described by a scattering 
matrix, or S-matrix, S(fc), relating incoming and outgoing amplitudes of propagat- 
ing plane waves of wave vector k = co/c (see Figure [lX| . Excellent pedagogical 
discussions, particularly for one-dimensional systems, can be found in.^^~^^ 



U Or 




Figure 1.1. Schematic description of the scattering matrix setup, showing the incoming (i^, ii) 
and outgoing (or, oi) amplitudes on the left and right of the scattering region. 

With obvious notations, the scattering S matrix is defined as: 

We consider the system to be invariant under time reversal, so that the matrix S is 
symmetric. Furthermore, it is unitary (S~^ = S"^)asa consequence of conservation 
of probability (for the Schrodinger equation), or of the intensity of the field (for the 
Helmholtz equation). This leads to the set of relations: 

|r|2 + = 1 (1.17) 

fr*+/f* =0. (1.19) 

These equations imply that det S = rr' — t^ = — f/f*. Since S is unitary, it can be 
diagonalized by a unitary transformation into the diagonal form: 



Defining the total phase shift, 5 (fc) = {(pi (fc) + ^2 (fc) ) / 2, we then have: 



(1.20) 



detS(fc) =e^^'^« = "F" ^^'^^^ 
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From the definition of the phase of the transmission amplitude, t = \t\e^^, and from 
( [1.21| , we obtain the relation S{k) = A:(fc) + n/2. A simple and elegant relation 
exists between the phase shift S{k) and the change of density of states p{k): 

p{k)-po{k) - i-Im^lndetS(/c) (1.22) 

where po{k) is the density of states for the free system, namely with zero potential 
for the Schrodinger equation in ( |1.3| , or = Cq for the Helmholtz equation ( |1.2| . 

The source of scattering in either Schrodinger or Helmholtz equations is the 
potential V{x) (or the varying refractive index n{x)). Should they vanish, the S- 
matrix reduces to the identity Now assume that they decrease fast enough so that 
we can enclose the scattering system (the "black box'' in Fig. inside a region 
of size L, much larger than the support of the scattering potential. Apply periodic 
boundary conditions, ip{0) = tp(L) and tp'(O) = ^'{L), at the boundary of the 
large box, noting that for large enough L, the physics is independent of the precise 
boundary conditions. For large enough L, 

ip\0) = xp'{L) ik{ii - oi) = ik{0Re^^^ - iRe-^^^) . (1.23) 

These algebraic relations may be written as a spectral condition: 

JkL (01 



deti^l-e^''^:^^jS{k)j =0. (1.24) 

Solving for S, leads to the following relation between the total phase shift defined 
in \1.21\ and the possible wave vectors: 

Noting that fc^^^ (L) = rcn/lL are the eigenmodes in the absence of scattering po- 



tential, namely for 3{k) = 0, we can rewrite (1.25> for two consecutive values of n 
under the form, 

(fc«+i-fcn)(L + ^) ^/r (1.26) 

Defining the density of states, or density of modes (DOM) as |0(/c) = 1/ (fc^+i — kn) 
leads to ( [L22| : 

PW-Mt) = i«^ (1.27) 



Recalling the definition in {1.9) of the counting function, M{k) = p{k')dk' , and 



its relation to the DOM in ( |1.1Q| ), we obtain: 

d{k) = 7tAN{k) = n {N{k) - No(fc)) (1.28) 

which relates the total phase shift of the S-matrix to a spectral quantity. Both \1.22\ 
and \1.2S\ are rather remarkable (and well-known) results since they express the 
fact that a measurement of the scattering data from a black box allows one to re- 
trieve its spectral information provided it is coupled to the external environment. 
Some further details are given in the next section. 
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1.4. Scattering matrix formalism | 7 
1 .4.1 . Derivation 2: Green's function and resoivent 

A more formal derivation, valid for all types of waves, is formulated in terms of 
Green's functions. A thorough description for the cases of interest here is pre- 
sented (among others) in Chapter 3 of/ Here we present some salient results as 
well as new expressions obtained using the Gelfand-Yaglom description. For prac- 
tical purpose, we consider the Schrodinger equation, but we keep in mind that 
these considerations apply equally to other wave equations. We denote by H the 
Hamiltonian of the scattering system and by Hq those of the free system, so that 
V = H — Hq defines the scattering potential. 

We define the Green's operators (resolvents) for the systems with and without 
the potential: 

G± = ^— , Gf= ^^/ ^ , £^0+ (1.29) 
E±ie — H ^ E±ie — Hq 



The relation to density of states p{E) defined in ( 1.1Q| follows from the identity: 



-lie 

Thus, we have 

p{E) = tr 3{E -H) = -^tr (G+ - G_) 

1 d , G-L 13,, G-L 

Here we have used the fundamental relation: 

tr In(operator) = In det(operator) (1.32) 

which is clear for diagonalizable matrices, but which also applies to the operators 
being considered here.^^'^^ In fact, the density of states should be defined relative 
to the free density of states: 

1 a f G G^°^ \ 



^ ^ In )^ (1.33) 



Ini dE 

det 



For a hermitean Hamiltonian H, the S-matrix is defined as^^ 
C G^°^ 

S = ^^ ^ S+ = S-i (unitarity) (1.34) 

G+ 



and so the result ( 1.22| follows 



These results are part of a vast literature on the subject generally known as the 
Krein-Birman-Schwinger formulation.^^"^^ It is based on the generalization of the 
resolvent operator G(z) = l/(z — H) to the complex z-plane. The determinant 
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of the S-matrix defined in ( |1.34| can be related to the resolvents defined in the 
presence (G) and in the absence (G^^^) of the scattering potential, 

IndetS(z) = lndet(z - H) - lndet(z - Hq) . (1.35) 

Thus, 

ImTr [G(£ + zO+) -Go(£ + /0+)] = (^-^^^^^^^i^)^ • (l-^^) 

Recalling that the density of states is given by Tcp{E) = — Im TrG(£ + zO+) leads 
again to ( |1.22[ ). This relation is a particular case of a more general set of relations 
known as the Krein-Birman-Schwinger relations which can be rewritten as, 

Tr[G(z)-Go(z)] = ___ IndetS(E) (1.37) 

or more generally, 

Tr[0(H)-0(Ho)] = lndetS(£) (1.38) 

where 0(H) is some regular function of the Hamiltonian H. When applied to 
= e~^^ , (1.38> gives the well known relations between the heat kernel F[t) = 
Tr and the S-matrix: 

/+00 
dEe-'^ [p(E) - poim (1-39) 
-00 

and between the zeta function ^h[^) = TrH~^ and the total phase shift: 

^„(s) - ^„„(s) = - l^yEE-'^ (1.40) 

We now introduce some terminologies used in different fields to describe re- 
lated quantities. From ( |1.34[ ), we have that 

detS(z) = det^^^ (1.41) 
z — Hq 

We can rewrite \135\ in the form In det S(z) = Tr ln(z — H) — Tr ln(z — Hq). Then, 
using the definition of the counting function A/'(£), we have 

iS{z) = J dJ\f{E) ln(z -^)- j d-^oi^) - £) (1-42) 

The quantity iS{z) = L(z) is sometimes called the Lyapunov function.^ Taking 
z = £ ± z£ with £ ^ 0+, we obtain 

S{E) = n{N'{E)-N'o{E)) (1.43) 

which states that the total phase shift measures the change of counting func- 
tion up to energy E in the presence of the scattering potential. The quantity 
A(z) = e^^(^) = det(z — H)/{z — Hq) is often called the Fredholm determinant 
or the Bottcher function in the mathematics literature. 

There is an interesting relation established by Herbert and Jones and by Thou- 
less,^^ which identifies the real part of the Lyapunov function 

7(£) = ReL(£ + ie) = J dN{E') In |£ - (1.44) 

as the spatial decay (or localization) length of the corresponding energy state. 
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1.4. Scattering matrix formalism | 9 

1 .4.2. Gelfand'Yaglom description and Transmission Probabiiity 

In fact, in addition to the density of states ^p{E), the Green's functions defined in 



the previous section 1.4.1 can also be used to compute the transmission probability 
T = \t\^. Such a combined description of spectral (density of states) and transport 
(total transmission probability) provides a very useful and powerful description of 
physical properties of quantum mesoscopic sy stems. ^'^^ The results can be sum- 



marized as follows. Define the ratios of Green's functions to free Green's functions 
as: W± = G± / G^^ . Then the transmission probability can be computed as the 
product of determinants: 

r(£) = |f(£)|^ = det(W+(£))det(W_(£)) = |detW+(£)|^ (1.45) 

while the density of states ( 1.33[ ) follows from the ratio of determinants: 



These relations are valid in any space dimension, but are sometimes difficult to im- 
plement in specific problems since they involve formal definitions of determinants. 
However, for one-dimensional scattering problems these calculations simplify dra- 
matically using the Gelfand-Yaglom theorem, which we now describe. 

Gelfand-Yaglom theorem: 

At first sight one might think that computing the determinant of a differential op- 
erator is a hopeless task, as there are infinitely many eigenvalues, which then need 
to be multiplied together. However, the Gelfand-Yaglom theorem is a remarkable 
result that allows one to compute the determinant of a (one-dimensional) differen- 
tial operator without computing any eigenvalues. ^^"^^ 

Consider the on e-dim ensional scattering problem with potential V{x), as illus- 
trated in the Figure 1.2. and the Hamiltonians H = — ^ + y, and Hq = — 
According to ([1.45[), the transmission probability is 



r(£) = |f(£)|^ = 



det (Ho - £ - it) 



(1.47) 



det {H-E- ie) 

The determinant can be computed in an efficient way using the GelTand-Yaglom 
theorem, which for our purposes can be stated as follows: 

Consider the Schrodinger eigenvalue problem Htp = \ip, on the interval x G 
— f / f / where L is, as previously, much larger than the range of the potential 
V[x). We apply Dirichlet boundary conditions, as this is relevant for the scattering 
problem, but the Gelfand-Yaglom theorem can be generalized to other boundary 
conditions.^^ Then the determinant of H can be computed as follows: solve the 
zero eigenvalue initial value problem'. Hu = 0, with u (^—j^ =0 and u^ ~ 
and similarly for HqUq = 0. Then the Gelfand-Yaglom result^^"^^ states that 
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Notice that this means that all one needs to do is to numerically integrate, with 
simple initial value boundary conditions, from one end of the interval to another, 
which is straightforward to implement numerically, as illustrated in figure pTzj For 
other boundary conditions there are analogus results. 



L 

^ = -2 




L 



Figure 1.2. The Gelfand-Yaglom theorem^^"^^ states that to compute the determinant of a 
differential operator (with Dirichlet b.c/s), one simply integrates numerically from one end 
of the interval to the other, with initial value conditions that the function vanishes and has 
initial slope 1. Then the value of the function u{L/2) at the other end of the interval is equal 
to the determinant. In fact, to be more precise this applies to ratios of determinants, as in 

Osl. 



To show how this is related to the transmission problem, we shift H by £ , as 
in (1.47>, and write E = k^. Then, near x = — L/2, where there is no potential, the 
solution of (H — fc^)z^ = 0, satisfying the initial conditions, is 



k ^ 2ik 



^^ik(x+L/2) 



-ik{x+L/2y 



(1.49) 



which has been decomposed into its ''right-mover'' and "left-mover" parts. Simi- 
larly, near x = +L/2, again the potential vanishes, so we can write 



^ ^z/c(x+L/2) ^ g ^-z/c(x+L/2) 



(1.50) 



for some coefficients A and B. To describe scattering from the right (for the prob- 
ability, the direction does not matter ^^'^^), we take fc ^ fc + ze, with e 0+, so 
that as L ^ 00 we select only the left-mover near x = —L/2. Near x = +L/2, 
this selects the left-mover, term, which is the incident wave. Thus, the 

transmission amplitude is 



r(fc) = 



2ik' 



-ikL/2 



1 



Be-ikL/2 

But Mo - 2^ (e^^'^+^Z^ 
Similarly, m(L/2) = Be''''^. So: 

m(L/2) 



2ikB 



(1.51) 



g-/J:(x+L/2)A gQthat Mo(L/2) 



^ e-*^ at X +L/2. 



2ik 



(1.52) 



from which the result \1A7) follows. 
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1.5. Transfer matrix formalism | 11 

1.5. Transfer matrix formalism 

In the previous section, we have advocated the description of wave propagation 
in complex systems using the S -matrix. It relates outgoing to incoming ampli- 
tudes in a unitary way, and provides a number of elegant expressions for spectral 
and transport quantities. However, in certain situations, such as for transport in 
one-dimensional problems, the transfer matrix formalism offers computational ad- 
vantages. For example, in layered systems it is useful to take advantage of the 
propagation of an incoming wave along a given fixed direction using the multi- 
plicative property of the transfer matrix. To see how this works, rewrite ( |1.16[ ) in 
an equivalent form, relating left-bound and right-bound amplitudes, instead of 
incoming and outgoing amplitudes (see figure [TX|: 



M (1.53) 

which defines the transfer matrix M. Note that det M = 1. For a symmetric poten- 
tial, namely for r' = r, M takes the interesting symmetric form: 

M=('^' -'/J\ ,1.54) 



{r/ty 1/f 

Remark: A great asset of the S-matrix and M-matrix formulations is that they 
are very general and they do not require to specify a wave equation or any lo- 
cal properties such as group velocity, dispersion at least until some more specific 
calculations. 

The main advantage of the transfer matrix description ( |1.54[ ) over the S-matrix 
is that it allows to perform the calculation of the total transfer M matrix for the 
case of a one-dimensional periodic, quasi-periodic (fractal) case since we have a 
multiplication structure. The transfer matrix can be viewed as a mapping trans- 
forming the wave after it pa sses through each scatterer or layer. Note first that 



the eigenvalues A± of (1.54> are solutions of — 2ARe(l/f) + 1 = 0. Since 
M is unimodular, A+A_ = 1 and propagating wave solutions correspond to 
—2 < 2Re {1/t) = TrM < 2, we can parameterize: 

TrM = 2cos^(fc) (1.55) 
which defines the so-called Bloch angle d(k). The eigenvalues of M(fc) are thus 
^±10 {k)^ There is an important relation between the Bloch angle and the total phase 
shift S{k) defined in \l.21 ). To establish it, we start from the relation 1 — e'^^^ 



It cos 6 between the two equivalent S and M formulations. Multiplying both sides 
by its complex conjugate, we obtain: 

s\n5{k) = ±|f|cos^(fc). (1.56) 
Using the expression ( |1.9[ ), we obtain a relation between the counting function and 
the Bloch angle: 

7tAJ\f{k) =S{k) = Arcsin {\t\ cos 0{k)) . (1.57) 
This expression is useful numerically. It allows one to obtain directly the counting 
function and the density of modes by derivation from a counting of the accumu- 
lated Bloch angle. Before entering into more details for the case of the Helmholtz 
equation, we present now a slightly different but yet closely related point of view. 
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1.6. Transfer matrices and a discrete Riccati equation 

For convenience, we consider in this section the case of the Schrodinger equation. 
Consider scattering above a single step, as indicated in the figure [T3] The wave- 
functions to the left and right of the step, located at x = Un, can be expressed as: 



n = (Ane'^-''- + Bne '^^""A , Xn<an 



where kn = ^E — Vn and fc^+i = ~ ^n+i- 



E 



(1.58) 



n+i = J (A„+ie''^"+i^"+i + B„+ie-''^"+i^"+i) , x„+i>fln (1.59) 



Vr 



n+1 



Vr 



n 



Figure 1.3. One segment of the discrete step potential for the discrete Riccati equation. 

Matching xp and xp' dX x = Un, we connect the coefficients of the left- and right- 
moving waves as: 

't;)KS9(B:;) 

where 

l^Jknkn+l '^V^nkfi+l 

For scattering from the left we set B^+i = and deduce the reflection and trans- 
mission amplitudes: 

'"=¥=-^ ' ^"=¥^=^ (1-62) 

Conservation of probability is satisfied at each step, since Ir^ p + | p = 1/ which is 
equivalent to | a:^ p — | jS^ p = 1, and which corresponds to the fact that the transfer 
matrix defined in \1.60\ has unit determinant. The transfer matrix equation ( |1.6Q[ ) 
can also be expressed as an evolution equation for the r^: 

rn+, = (1.63) 

which we recognize as a discrete form of a Riccati equation. 
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1 .6.1 . Continuous limit of discrete Riccati equation 

In order to take the continuous limit of a smoothly varying potential, for scattering 
above the potential, we think of the potential as broken into a series of discrete 
steps and successively apply the transfer matrix formalism of the previous section 
to pass from one step to the next. It proves convenient to use a slightly modified 
phase convention than in ( |1.59[ ): in the Xn segment, a^^i < Xn < Un we define 

(1.64) 

and in the segment, Un < x^+i < a^-^-i we define 

(1.65) 



This phase convention refers the phase to the left-hand end of the segment, and 
keeps a total accumulated phase from previous segments. Then the transfer matrix 
relation {i.60\ becomes: 



Dn+1 J \Pn ^nj J 

where now dcn is real: 

^ ^n+l+}n_ ^ ^ fcn+^^-2.T»*:A-M) (1.67) 

Clearly we still have unitarity: — li^np = 1/ for all n . Def ining a reflection 
amplitude at step n as = / A^, the evolution equation ( 1.63 > can be written as 



del ^ K K _ Pn 
^n+l —^n — : 1 — -z- 

^ lyt 



n* 'n 



(1.68) 



where we have used the fact that the ocn are real with this phase convention. 

To take the continuous limit, we take the Un to be equally spaced by distance A, 
and note that 

k„+t-kn 11 Vn+l-Vn ^ lV'{x) 

A ^ 2A ^ lQ{x) ^ 

1 1 

—— (1.70) 



kn+1 + k„ 2Q{x 
where Q(x) — a/E — V{x). Therefore, 

kn+l-k„ \ 2il«kj{aj-aj_,) ^ Q' jx) ^2irQix')dx' q y-^. 

A a* A \kn+i+knj 2Q{x) 

and similarly for the conjugate. Furthermore, ^n/^n — C'(A) — > 0. Therefore, the 
continuum limit of the discrete Riccati equation l |1.68| is 

f'(x) = g (e^'TQ- f2(x)e-2^/'Q) (1.72) 
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This is straightforward to implement numerically: integrate for scattering from 
the left with the initial condition f(— oo) = 0, yielding the reflection amplitude as 
f(+oo). 

For comparison we briefly recall the derivation of this continuous Riccati equa- 



tion (1.72) for one-dimensional scattering/^ Re-write the wavefunction ip{x) in 



terms of WKB adiabatic solutions with coefficients A{x) and B{x) as: 

= ^^^=y (A(xyre + B(x)e-rQ) (1.73) 

rp'ix) = -^gl= {Aix)e'rQ-Bix)e-'rQ) (1.74) 



where, as above, Q{x) = ~ V{x). The decomposition (1.74> requires the fol- 
lowing relation between the (complex) coefficient functions: 

A'{x)\ _ Q'{x) ( e-^irQ\ (A{x) , ^^^^^ 



B'{x)J 2Qix) [e^ifQ J \B{x) 

Then defining the reflection coefficient amplitude f{x) = ^|^/ we see that it satis- 
fies the differential equation 

in agreement with the limit \1.72) of the discrete transfer matrix approach. 
Comments: 

(1) We should distinguish the Riccati form \1.76\ from another, perhaps more fa- 
miliar, realization of the Schrodinger equation as a Riccati equation: define a 
''momentum'' function in terms of the derivative of the wavefunction xp{x) as 
P ~ ~^Ix Schrodinger equation for tp{x) becomes a Riccati equa- 
tion for p{x) 

= iQ-i/ (1.77) 

A simple computation shows that the two Riccati equations \l.76\ and \l.77\ 
are completely equivalent, with the functions and f(x) related by: 

(2) It is worth noting that when integrating a Riccati differential equation, such as 
y'{x) = — c(x) y^(x), it often proves numerically useful to use a Mobius 
integrator^^ for step A 

rather than a conventional Numerov integrator. Clearly, the discrete Riccati 
equation ^1.63) is precisely of this form, when we make the phase choice ( 1.67[ ). 
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1.7. Transfer matrix description for the Helmholtz equation 

The purpose of this section is to implement the previous methods in the case of the 
propagation of electromagnetic waves in one-dimensional (ID) waveguide struc- 
tures built out of complex media. We first reiterate the basic equations and then 
consider specific cases of layered dielectric media : free space, Fabry-Perot struc- 
ture, periodic structures (photonic crystals), random systems and single impurity 
Finally, we discuss in more detail the case of a Fibonacci potential. In all these 
cases, we obtain spectral quantities (transmission, counting function and density 
of modes) and the steady state (stationary) local structure of the electric field. 

The use of the transfer matrix method for solving Maxwell equations in pho- 
tonic systems essentially converts them into a set of finite difference equations in 
real space and then rearranges those equations into the form of a transfer matrix.^^ 
In ID structures, transfer matrices relate the electric and magnetic fields in one 
plane those in an adjacent plane. The method can serve as an approximation for 
any dielectric stratified medium with e = e (z) ,}i = }i (z). In the case of an arbi- 
trary finite structure of piecewise constant potential, the transfer matrix method is 
an exact solution. 

The general solution for the transverse electric (TE) electric field Ex propagating 
along the z-direction in a ID structure is obtained from the solution of the wave 
equation c^d^Ex — d^Ex = 0. It takes the form 

Ex (z, t) = (zo) + Ei^^ (zo) e^'(-^z(^)^-a^O (1.80) 

where E^^^ (zq) and (zq) denote rightbound (R) and leftbound (L) waves 
steady-state amplitudes evaluated at the origin zq of the structure. Setting fc(z) = 
n (z) fc, we can write for a monochromatic wave of frequency co, 

(z) - Ef ) (zo) e'^^'> + E^P (zq) . (1.81) 

Using the definition \153) of the transfer matrix M, we obtain 

where Zq and z^^^ are the coordinates of the two outer boundaries of the structure. 
In non-absorbing dielectric structures considered here, conservation of energy im- 
plies that M is a unitary 2x2 complex valued matrix of unity determinant. In a 
layered structure Zq < Zi < • • • < z^^^, and denoting M^.^z-^^ the transfer matrix 
associated to a layer, we have a general multiplication scheme: 

^z,^z,,, = M,^^,_^^,^^, • . . M,,^,, M,,^,, . (1.83) 

For practical purposes we consider from now on, the scheme represented in Fig- 
ures 



1.4. 



and 1.5. is a bit more intricate. Namely, each M^.^z.^^ appears as a 
product of sub-matrices responsible either for an interface boundary between two 
adjacent dielectrics or for the free propagation through a homogeneous dielectric 
slab. 

We will now restrict ourselves to the more specific case of a binary layered 
structure built out of N slabs of two types of dielectrics A and B, where g are 
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the respective refractive indices in type A and B slabs of respective corresponding 



thicknesses ^^,5- The notations for this case are given in figures [lA] and 1.5. The 
solution for the electric field in a single slab of type A or B can be written as: 



are t 

at the left boundary of slab A, B 



(1.84) 



where and are the right-bound and left-bound electric fields evaluated 



Remark: As stated in the Remark in section 1.5 the M (or S)-matrix formalism 
does not require precise knowledge of a local wave equation inside the structure 
and can be formulated without the input of a local group velocity (even ill-defined 
in the present case due to the discontinuity of the refractive index in the layered 
structure). But at some point the writing of relations \1.S1\ and \1.S4\ rely on the 
definition of a wave vector and a group velocity : here we assume free propagation 
in each individual slab. 




A 










Figure 1.4. Setup of notations for the continuity relations at an interface ABA (upper fig- 
ure) and for the free propagation through a slab (lower figure). Upper part: fields at the 
interfaces AB and BA are noted. Lower part: fields at the right and left ends of slab A and B 
are noted. 



An interface 1 = A ^ B between adjacent slabs A and B is defined in figure 



1.4. The boundary conditions at a dielectric interface are imposed through conti- 
nuity relations: 



(L)| 
A M+, 



— Ct, -|- Er II, 



(continuity) 



(1.85) 



and 



Ua (sjf ^ - = riB (Ef^ - (derivative-continuity) (1.86) 

The phase accumulation which accounts for free propagation through a slab A or 
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Figure 1.5. Setup of notation for the ID transfer matrix model for a binary layered struc- 
ture. Here z is the ID coordinate, dj^^B is the slab thickness, n^^ g is the refractive index, and 
E^'^ is the right-bound and left-bound electric field amplitudes. 



B is: 



(R) 



,(L) 



(R) 



,(L) 



BdA,B^ = f(^) p-i^A,B 

where we have defined, 

5 a fiAd-Ak , Sb fiBdBk. 
It is convenient to set the layer thickness so that n^d^ = ngdg which leads to 
^A = ^B- With the change of variables: 

:(+) = ^ F^) 

^A,B 



(1.87) 



(1.88) 



r(-) 
-A,B 



^A,B 



\^A,B ^A,B ) ' 

continuity equations ( |1.85[ ) and ( |1.86[ ) can be rewritten as: 



(1.89) 



-A,B 



1+,- 



= E 



( + ) 



^BA 



(1.90) 



This leads to the definition of four transfer matrices defined by the two sets of 
relations. 



(1.91) 





and 




A^ 



B^ 



(1.92) 



With the aid of equations ( 1.88 > and (1.89 >, we obtain four real valued matrices :^^'^^ 



Tab — 



1 



BA 



Tba = 




= T 



AB 



(1.93) 
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and 



^ ^^cos5^,,-sin5^,, 
\sin5A,B cos5a,b ' 



so that the total transfer matrix Mjv of N dielectric slabs is 

N 



To illustrate the efficiency of this systematic approach to the calculation of Mjv as 
a product of an appropriate sequence of the four sub-matrices ( |1.94[ ), we consider 
the matrix Ms describing an arbitrary structure inserted between slabs of type A 
(to preserve unitarity) as: 

Ms{[A][ABABAABA][A]) = TaTabTbTbaTaTaTabTbTbaTaTabTbTba (1-96) 

Once the transfer matrix is obtained, an algebraic transformation between ( |1.16[ ) 
and \1.53\ allows to derive the scattering matrix 

(|;)-.(f) 

which relates incoming to outgoing fields. From these expressions, we deduce the 
transmission and reflection coefficients T {co) and Tl{a)), the density of modes 
p [co) and the counting function M {co). 



Transmission and reflection spectra : The transmission and reflection spectra 
are calculated for the case of an incident wave from the left (see |1.98[ ) and spectra 
are normalized by the incoming plane wave electric field amplitude, namely. 



E'^y = t{cv) 



4'^ = r{co) 
£(^) - 




(1.98) 



Using the transfer matrix Mtot, defined by 



M, 



tot 



i_y 1 , yields: 



r{cv) 
nicv) ■ 



4|M(a;)|^ 4 
||M(a;)||2+2 ~ ||M(a;)||2+2 



||M(a;)f-2|M(a;)| ||M(a;)f 



||M(a;)f+2|M(a;)| ||M(a;)f+2 

where ||Mf = Mf^ + + M^^ + M^^, and M = [M/y]. 



Ticv) 



(1.99) 
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Density of modes spectrum calculation: The density of modes (DOM), p (a;), 
can now be calculated from either the scattering or transfer matrix, utilizing the 
total phase shift S{co) or the phase of the transmitted amplitude oc (co) defined in 
section ( |1.4| . For a ID system of length L, the normalized DOM given by \1.22) is 
p {co) = ^-^^ = i^^^- The normalization with the length L is convenient as 
it keeps the DOM of a free system equal to one, and simplifies the evaluation of 
the relative enhancement of p (to) for a structured system compared to a free one. 
This normalization also allows, in the case of aperiodic inhomogeneous systems, 
to give a meaning to otherwise ill defined quantities such as the effective group 
velocity and wave vector. The calculation oi p (co) using either scattering \1.22\ or 
transfer^^ matrices is thus based on, 

(1.100) 



Amplitude of the local electric field : The electric field amplitude and inten- 
sity CO — z maps are obtained by calculating E (to) for each z coordinate inside the 
structure, using two transfer matrices: the total transfer matrix Mtotr and the par- 
tial transfer matrix from the left end to the point of interest My. Firstly, r and t are 



calculated as before using Mtot/ 
intensity in slab ; using: 



and then we can use My to calculate the field and 




(1.101) 



The net electric amplitude E{z)^^ at a point z, is the coherent amplitude sum of 
multiply reflected right-bound and left-bound waves, incorporating interference 
effects inside the system. This process is repeated for all ; in order to create the 
full electric field E^^^ (a;,z) and intensity distribution I^^^ {co,z) maps. The spatial 
resolution can be enhanced by numerical subdivisions of the slabs. This calculation 
can be performed for the cases of an incident wave either from left, right or from 
both sides. 



1.8. Illustrative Examples of Layered Systems 
1.8.1. Free space 

In this illustration, the entire structure is composed of type A slabs. The transmis- 
sion spectrum and the normalized density of modes are shown in figure \1.6.\ . As 
expected, both quantities are equal to one and independent of co. The electric field 
amplitude map is displayed in figure ( |1.7.[ ). The map shows the fundamental si- 
nusoidal oscillations of the traveling plane wave, as seen in the cross section of the 
map on the right part of the figure. 
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Figure 1.6. The normalized density of modes (left part) and the transmission spectrum 
(right part), for free space (no scattering structure). 




1 .8.2. Fabry-Perot structure 



The all dielectric Fabry-Perot resonator is realized here by placing two type B slabs 
separated by several type A slabs, and for illustration we have chosen the refractive 
index to be three times larger than n^. The transmission spectrum and the 
normalized density of modes are presented in figure dLSl. 
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Figure 1.8. The normalized density of modes (left part) and the transmission spectrum 
(right part), for a Fabry-Perot resonator. 



Although this simple realization serves as a poor resonator, transmission peaks 
and dips are distributed according to the Fabry-Perot resonance condition = 
{^}n=oi2 ' where the density of modes is enhanced above the free one at reso- 
nances and attenuated off-resonance. This is, of course, due to interference effects 
inside the resonator as displayed by the corresponding electric amplitude map in 
figure \1.9.) . 

The electric field amplitude map is identical to the free space map with the 
additional "selection rules'' for the Fabry-Perot cavity (''allowed'' and "forbidden" 
modes). The two different cross sections of the map shown in the right part of 
figure \i.9.) make clear the strong destructive interference effects taking place in 
'off resonance' frequencies inside the structure. 




Figure 1.9. Left panel : the electric field amplitude map for a Fabry-Perot resonator. Right 
panel: a cross-section of the map along the entire stack at a frequency of an allowed mode 
(2L = 6A ) in blue and of a forbidden mode (2L = 7.5A ) in red. 



1 .8.3. Periodic structure - phiotonic crystai 

The transfer matrix formalism introduced in section ( |1.5[ ) is particularly convenient 
to study the case of a periodic layered system.^^"^^ In this case, the calculation of 
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the total transfer matrix {Mceii{co))^ for a sequence of N slabs, simplifies consid- 
erably thanks to the Cayley-Hamilton theorem which states that for unimodular 
matrices, 

- (TrM) M + 1 = (1.102) 
This allows to easily obtain in a recursive way, powers of the transfer matrix 
M^(fc) = L7^_i(cos^)M(fc) - LI„_2(cos^) (1.103) 

where Un{cosO) = ^^^ginV^^ Chebychev polynomials of the second kind.^^ 

A direct calculation shows that 

TrM^ = cos ne{k) = T„(cos^) (1.104) 

where Tn (cos 6) are the Chebychev polynomials of the first kind. It is also straight- 
forward to obtain the expression of the transmission coefficient Tjv = |f]vp for a 
Bragg system made of N identical layers each described by the same transfer ma- 
trix M. From ( 1.103 >, we have the expression of M^, which by definition is of the 



form ( |1.54[ ), namely. 



1 



Un- 



t 



(1.105) 









Chebyshev 
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Figure 1.10. Periodic Bragg system. The band structure is calculated using both the Cheby- 
chev polynomials jl.l06^ and via the transfer matrix approach ( [1.100^ . Right part - a close 
up on the long wavelength band edge mode. 



1 


2 








tN 




tN 



Since det = 1, then 



= 1, so that 
1 



1 I I r |2 / sin NO \ 



(1.106) 



This is a straightforward way to numerically obtain the well known Bragg spec- 
trum (see Figs. 1.10.| and |l.ll.[ ). 
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Figure 1.12. The transmission spectrum (top part) and the normalized density of modes 
(bottom part) for a 40 slab photonic crystal. The figures on the right shov\^ a close-up on the 
long v^^avelength band edge mode 



Remark 

The transmission coefficient Tj^{0) for a finite periodic system (10 unit cells in 
length) were calculated using both a direct transfer matrix calculation and the re- 
lation ( 1.22[ ) and more directly ( 1.106[ ). Both are plotted in figure ( [1.10. ) and the 



agreement is perfect. From the relation ( |1.106[ ), it appears clearly that the reso- 
nances are given by the condition = 0. Note that this is specific to the peri- 
odic case and cannot be generalized to other more complex structures, for which 
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the general relation must be used to obtain the spectrum. 

The dielectric photonic crystal analyzed from here forth is composed of alter- 
nating type A and B slabs with 20 repetitions of the unit cell [AB]. The refractive 
index is chosen to be = |n^. The transmission spectrum for this 40 slab pho- 



tonic crystal is shown in Figure 1.11. The periodicity gives rise to transmission 
bands and photonic band gaps centered at the ^ condition, as expected. A closer 
look at the transmission and the density of modes spectrum near the first stop band 
is depicted in Figure [LIZ 



The sharp band edge modes with a 50-fold enhancement in the DOM compared 
to a free system can be seen in the close up figures. The number of resonant modes 
corresponds to the number of unit cells. 

Remark : All spectral quantities presented here were calculated using the total 
transfer or the total scattering matrix, but are are completely equivalent to ( |1.106[ ). 

The electric field amplitude map corresponding to this photonic crystal is dis- 
played in Figure 1.13.[ as well as a close up look at the band edge modes. Figure 




0.32 0.33 

6/71 



0.67 0.68 
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Figure 1.13. The electric field amplitude co — z map for a 40 slab photonic crystal. Spatial 
boundaries of each slab are indicated by dashed gray lines. Bottom figures are a close-up 
on the long (bottom left figure) and short (bottom right figure) band edge modes. 



1.13. 1 shows that the band edge modes also have an enhanced electric field distri- 
bution peaks - up to 12 times that of the incoming field, due to constructive inter- 
ference effects and transient energy buildup. It is also visible that while the electric 
field ''hot-spots'' in the long wavelength band edge mode tends to be located in 
one type of slab, it tends to be located in the other type for the short wavelength 
band edge mode. This feature is presented more clearly in Figure p^A] 
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Figure 1.14. the electric field amplitude cross section for the long and short wavelength 
band edge modes in a 40 slabs photonic crystal. Electric field of the short and long wave- 
length modes concentrate on different slab types. 



1.8.4. Random structure 



Disordered layered structures have been studied in very great details as a success- 
ful and easy to implement model for ID realization of Anderson localization^"^ 
(and references therein) and other related complex systems. It w^ould be a hope- 
less task to list all relevant works here and it is also not our purpose to study this 
class of systems in greater details. We discuss them in this review^ mostly for the 
sake of completeness and also to present results for the electric field amplitude 
map not so often discussed in the existing literature. 

A dielectric random system is generated here by a random sequence of type 
A and B slabs. Weak disorder is obtained by setting the refractive index tlb to be 
10% larger than n^, and for strong disorder, the refractive index contrast is set 
as riB = 2n^. The transmission spectrum and the normalized density of modes 
for a w^eakly disordered structure is displayed in figure |1.15. It can be seen in 
the right part of figure [1.15. that a weak disorder can attenuate the transmission 
spectrum significantly, but as seen in the left part of figure |1.15. the density of 
modes shows no localized modes - the localization length in this case is greater 
than the total length of the structure, so the modes appear "smeared'' into each 
other. The transmission spectrum and the normalized density of modes for strong 
disorder is shown in figure [lT6] 
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Figure 1.15. The normalized density of modes (left part) and the transmission spectrum 
(right part), for a weakly disordered structure 
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Figure 1.16. The normalized density of modes (left part) and the transmission spectrum 
(right part), for a strongly disordered structure 



In contrast to the previous case of weak disorder, the transmission for strong 
disorder is inhibited for all frequencies except for a set of discretely distributed 
modes. The density of modes shows the same discrete resonant frequency pattern 
together with a significant enhancement of the DOM campared to a free system. 
The electric field amplitude map for weakly and strongly disordered structures is 
shown in figure pTTzj 

The modes of a weakly disordered structure shown in subplot (a) of figure 



1.17. are smeared spatially and the spatial distribution of electric field extends to 



the edges of the stack. For the case of strong disorder in subplot (b) of figure pTTz!} 
the electric field amplitude map displays vanishingly small values for all frequen- 
cies except some discrete isolated modes as the one shown in subplot (d) of figure 



1.17. These modes are spatially localized with a localization length smaller then 



the structure length. The electric field amplitude is greatly enhanced in the reso- 
nant modes - about 70 times that of the incoming wave - also due to constructive 
interference and transient energy buildup. 
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1 .8.5. Fibonacci structure 

Among the large variety of complex structures, quasi-periodic distributions play 
a special role. One of the many reasons for this ongoing and continuous interest, 
is that quasi-periodic structures are deterministic (non random) but nevertheless 
display interesting qualitative similarities (but quantitative distinct features) with 
random structures, e.g. localization of the modes, insulating behavior, structured 
spectra,^'^ 

Moreover, quasi-periodic ID structures can be viewed as realizations of quasi- 
crystalline order and there are interesting candidates to study transport and thus 
understand the insulating behavior of these structures. In this short review, we 
consider a two-word Fibonacci and a Cantor set structure as working examples to 
test our approach to obtain spectral properties using the S-matrix as in ( 1.22[ ). An 



interesting point about the Fibonacci case is that while this real space structure is 
not fractal (it does not display a discrete scaling symmetry leading to a self -similar 
order), the corresponding spectral quantities display a self-similarity property yet 
to be fully defined. ^'^'^O'^^-^^ 

The literature on quasi-periodic distributions is extremely rich. Here we 
present some references which may be considered as leads to a more complete 
bibliography. The initial interest in such problems traces back to works by Dyson 
and Schmidt.^ For a recent review on spectral properties see,^ and.^ Applications 
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of quasi-periodic structures to electronic transport properties and to electromag- 
netic wave propagation are numerous. More closely related to the topics covered 
in this book are.^'^^'^^'^^'^^ 

The Fibonacci structure can be defined in many different ways. Generally, con- 
sider an alphabet A= { A, B, C, • • • } composed of different letters, and a substitu- 
tion rule f {A), so that 

C{ABC---)=aA)maC) ••• (1.107) 

Then, a Fibonacci sequence is Fn = ^^{Wq), where Wq is some initial word. Her, 
we consider an alphabet composed of two letters A and B namely slabs of type A 
and type B arranged according to a Fibonacci sequence f (A) = BA, f (B) = A and 
f^ = f^(B), namely, 

A^BA^ ABA BAABA ABABAABA ^ • • • (1.108) 

Alternatively, the same structure can be defined by the first two words, and a build- 
ing rule stating that the third word and on is composed of the sum of the previous 
two words, namely, 

Fo = B; fi = A; Fj^^ - [Fj-2Fj-i] (1.109) 



An illustration of the resultant structures is depicted in figure 1.18.[ Setting the 
refractive indices to be = the transmission spectrum and the normalized 
density of modes for 10th generation Fibonacci structure are shown in figure [lT9^ 




A B A B A A B A^H B A ABAABABAABA 



Figure 1.18. A schematic view of generations 0-6 of the Fibonacci structure used in the 
calculations 

The spectrum shows two distinct band gaps, and some spectrally isolated band 
edge modes with enhanced DOM, arranged according to a self-similar "fractal 
like'' structure. Note that the DOM does not obey an exact self-similarity, but only 
an approximate one. However, it has been established^'^^ that the energy spectrum 
of an infinitely long Fibonacci structure is multi-fractal from the cantor set family 
(see figure [L23]). To explore this self similarity, we consider the transmission spec- 
trum for two different generations of the structure. The spectrum is depicted in 
figures [Tm] and [L2L| 
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Figure 1.19. The normalized density of modes (left part) and the transmission spectrum 
(right part), for the 10th generation Fibonacci structure. The structure of the stack is illus- 
trated above the two spectra. 




Figure 1.20. Transmission spectrum and stack structure for 8th generation (left part), and 
14th generation Fibonacci structure (right part). Structures are illustrated above the spectra. 
The self similarity of the spectra is evident in the close-up on the middle part of the two 
spectra. 




5/71 5/71 5/71 



Figure 1.21. the counting function J\f (cu) for a 10th generation Fibonacci structure. The 
middle and right figures are a gradual zoom on the counting function depicted on the left 
figure. A qualitative self similarity is observed. 



The similarity between the middle section of the transmission spectra for dif- 
ferent Fibonacci structure can be seen in the close up figures (the middle spectral 
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feature is the same for very different band widths). This result initially obtained 
jj^30,3i reproduced here. The electric intensity map for the 8th generation Fi- 
bonacci structure, and a close up look at some interesting modes is shown in figure 
1.22. The plots in figure \L22] show that the electric field intensity is vanishingly 
small for almost all frequencies except for some discrete modes. For some of these 
modes the electric intensity spatial distribution is localized and (surprisingly) sym- 
metric around the middle of the stack, and self similar as it resembles the Cantor 
set structure - see figure [1.23. The reason for the symmetry is believed to be the 
palindromic nature of the Fibonacci sequence. This somewhat odd electric field 
spatial structure is probed due to the enhanced intra-slab resolution of the calcula- 
tion. Here too there is a strong enhancement of the electric field amplitude in the 
band edge modes - up to 60 times that of the incoming field. 
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Figure 1.22. The electric field intensity map for the 8th generation Fibonacci structure. 
Left- full spectrum. Other plots: closeup zoom on some isolated modes. Spatial localization 
and spatial self-similarity are apparent. 
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Figure 1.23. A schematic view of generations 0-6 of a triadic cantor set. 



1 .8.6. Fractal Cantor set structure 

The Cantor set structure is realized here by arranging type A and type B slabs ac- 
cording to the triadic Cantor sequence (see figure [1.23. [ ), and setting the refractive 
index to be 50% larger than n^. The transmission spectrum and the normal- 
ized density of modes for a 4th generation Cantor set structure is shown in figure 
1.24. 1 Figure [1.24. shows that as for the random and the Fibonacci structures, this 
structure too possesses a jagged energy spectrum with discrete resonant modes ac- 
companied by a DOM enhancement. A close up look at the transmission curve of 
one of these modes, and the electric field intensity map for the same mode is given 
in figure 1.25. 
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Figure 1.24. The normalized DOM (left part) and the transmission spectrum (right part), 
for the 4th generation triadic Cantor set structure. The structure of the stack is illustrated 
above the two spectra. 
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Figure 1.25. A close-up on the transmission curve (left) and electric field intensity map 
(right) for a resonant mode in a 4th generation Cantor set structure. 



Figure 1.25. shows that as for the Fibonacci structure, the resonant modes in a 
Cantor set structure are isolated spectrally, and localized spatially. How^ever, the 
spatial localization in the Cantor case is geometric in nature - the electric field spa- 
tial distribution is enclosed in the middle third of the structure, while the smaller 
Cantor motifs at both sides serve as dielectric mirrors. 



1 .8.7. Defect modes in photonic crystai structure: 

The Fibonacci and Cantor set structures are examples for quasi-periodic structures 
which give rise to deterministic modes isolated spectrally and localized spatially. 
The "work horse'' of this kind of modes is the defect mode - a defect (say an [AA] 
segment) is inserted to the center of a photonic crystal to produce a mid gap mode, 
localized spatially inside the defect.^ The transmission spectrum and the normal- 
ized density of modes for a photonic crystal hosting a single impurity is shown in 
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figure 1.26. 
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Figure 1.26. The normalized density of modes (left part) and the transmission 
(right part), for a photonic crystal hosting a single defect. 
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Figure 1.27. A close-up on the transmission curve (left part) and the electric field intensity 
map (right part) for the resonant defect mode in the photonic crystal structure 



A close up look at the defect mode transmission and its electric field amplitude 



map, is depicted in figure 1.27. The defect mode is indeed isolated spectrally and 
localized spatially. How^ever, if unintentional defects also exist, additional defect 
modes can exist and may "smear'' each other spatially and spectrally. 



1.9. Discussion 

This chapter has presented the basic theoretical tools needed to describe the spec- 
tral and transport properties of complex layered systems. These techniques are 
quite standard in the fields of electrodynamics and quantum mechanics, but in- 
teresting new effects emerge when applied to transport in complex layered struc- 
tures. The explicit numerical results for various types of layered samples show^ 
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that deterministic a-periodic structures can serve as a deterministic complex cav- 
ity, hosting localized modes which are also narrow-band and spectrally isolated. 
In the Fibonacci dielectric structure, these modes also appear to be spatially self 
similar. This suggests promising applications for the physics of spontaneous and 
stimulated emission, and these initial results deserve to be explored further, exper- 
imentally and numerically. 

This work was supported by US DOE under grant DE-FG02-92ER40716, and by 
the Israel Science Foundation Grant No.924/09. We thank E. Gurevich for helpful 
discussions and for relevant remarks. 
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